This notebook demonstrates deposition of an SDS adsorption layer on a non-spherical AFM tip model.
%load_ext autoreload
%autoreload 3
import ase.io # here used for reading pdb files
from ase.visualize import view
from ase.visualize.plot import plot_atoms # has nasty offset issues
from cycler import cycler # here used for cycling through colors in plots
import datetime
import fabric # for pythonic ssh connections
from fireworks import LaunchPad, Firework, Tracker, Workflow
from fireworks import FileTransferTask, PyTask, ScriptTask
# FireWorks functionality
from fireworks import Firework, LaunchPad, ScriptTask, Workflow
from fireworks.user_objects.firetasks.templatewriter_task import TemplateWriterTask
from fireworks.user_objects.firetasks.filepad_tasks import AddFilesTask, GetFilesTask, GetFilesByQueryTask
from jlhfw.fireworks.user_objects.firetasks.cmd_tasks import CmdTask
from fireworks.utilities.filepad import FilePad # direct FilePad access, similar to the familiar LaunchPad
import glob
import gc # manually clean up memory with gc.collect()
import gromacs # GromacsWrapper, here used for evoking gmc commands, reading and writing .ndx files
# from io import StringIO, TextIOWrapper
import io
import itertools # for products of iterables
import json # generic serialization of lists and dicts
import jinja2 # here used for filling packmol input script template
import jinja2.meta # for gathering variables in a jinja2 template
from jlhfw.fireworks.user_objects.firetasks.cmd_tasks import CmdTask # custom Fireworks additions
import logging
import matplotlib.pyplot as plt
import MDAnalysis as mda # here used for reading and analyzing gromacs trajectories
import MDAnalysis.analysis.rdf as mda_rdf
import MDAnalysis.analysis.rms as mda_rms
from mpl_toolkits.mplot3d import Axes3D # here used for 3d point cloud scatter plot
import miniball # finds minimum bounding sphere of a point set
import nglview
import numpy as np
import os, os.path
import pandas as pd
import panedr # reads GROMACS edr into pandas df, requires pandas and pbr
import parmed as pmd # has quite a few advantages over ASE when it comes to parsing pdb
from pprint import pprint
import pymongo # for sorting in queries
import scipy.constants as sc
import subprocess # used for evoking external packmol
import sys
import tempfile
import yaml
GromacsWrapper might need a file ~/.gromacswrapper.cfg with content
[Gromacs]
tools = gmx gmx_d
# gmx_mpi_d gmx_mpi_d
# name of the logfile that is written to the current directory
logfilename = gromacs.log
# loglevels (see Python's logging module for details)
# ERROR only fatal errors
# WARN only warnings
# INFO interesting messages
# DEBUG everything
# console messages written to screen
loglevel_console = INFO
# file messages written to logfilename
loglevel_file = DEBUG
in order to know the GROMACS executables it is allowed to use. Otherwise,
calls to gmx_mpi or gmx_mpi_d without MPI wrapper might lead to MPI
warnings in output that cause GromacsWrapper to fail.
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger()
logger.setLevel(logging.INFO)
ParmEd needs to know the GROMACS topology folder, usually get this from
envionment variable GMXLIB:
def find_undeclared_variables(infile):
"""identify all variables evaluated in a jinja 2 template file"""
env = jinja2.Environment()
with open(infile) as template_file:
parsed = env.parse(template_file.read())
undefined = jinja2.meta.find_undeclared_variables(parsed)
return undefined
def plot_side_views_with_spheres(atoms, cc, R, figsize=(12,4), fig=None, ax=None):
"""
Plots xy, yz and zx projections of atoms and sphere(s)
Parameters
----------
atoms: ase.atoms
cc: (N,3) ndarray
centers of spheres
R: (N,) ndarray
radii of spheres
figsize: 2-tuple, default (12,4)
fig: matplotlib.figure, default None
ax: list of three matploblib.axes objects
"""
logger = logging.getLogger(__name__)
atom_radii = 0.5
cc = np.array(cc,ndmin=2)
logger.info("C({}) = {}".format(cc.shape,cc))
R = np.array(R,ndmin=1)
logger.info("R({}) = {}".format(R.shape,R))
xmin = atoms.get_positions().min(axis=0)
xmax = atoms.get_positions().max(axis=0)
logger.info("xmin({}) = {}".format(xmin.shape,xmin))
logger.info("xmax({}) = {}".format(xmax.shape,xmax))
### necessary due to ASE-internal atom position computations
# see https://gitlab.com/ase/ase/blob/master/ase/io/utils.py#L69-82
X1 = xmin - atom_radii
X2 = xmax + atom_radii
M = (X1 + X2) / 2
S = 1.05 * (X2 - X1)
scale = 1
internal_offset = [ np.array(
[scale * np.roll(M,i)[0] - scale * np.roll(S,i)[0] / 2,
scale * np.roll(M,i)[1] - scale * np.roll(S,i)[1] / 2]) for i in range(3) ]
###
atom_permut = [ atoms.copy() for i in range(3) ]
for i, a in enumerate(atom_permut):
a.set_positions( np.roll(a.get_positions(),i,axis=1) )
rot = ['0x,0y,0z']*3#,('90z,90x'),('90x,90y,0z')]
label = [ np.roll(np.array(['x','y','z'],dtype=str),i)[0:2] for i in range(3) ]
# dim: sphere, view, coord
center = np.array([
[ np.roll(C,i)[0:2] - internal_offset[i] for i in range(3) ] for C in cc ])
logger.info("projected cc({}) = {}".format(center.shape,center))
color_cycle = cycler(color=[
'tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple'])
circle = [ [ plt.Circle( c , r, fill=False, **col) for c in C ] for C,r,col in zip(center,R,color_cycle) ]
margin = 1.1
# dim: view, coord, minmax (i.e., 3,2,2)
plot_bb = np.rollaxis( np.array(
[ np.min(center - margin*(np.ones(center.T.shape)*R).T,axis=0),
np.max(center + margin*(np.ones(center.T.shape)*R).T,axis=0) ] ).T, 1, 0)
#plot_bb = np.array( [ [
# [ [np.min(c[0]-margin*R[0]), np.max(c[0]+margin*R[0])],
# [np.min(c[1]-margin*R[0]), np.max(c[1]+margin*R[0])] ] ) for c in C ] for C,r in zip(center,R) ] )
logger.info("projected bb({}) = {}".format(plot_bb.shape,plot_bb))
if ax is None:
fig, ax = plt.subplots(1,3,figsize=figsize)
(ax_xy, ax_xz, ax_yz) = ax[:]
logger.info("iterators len(atom_permut={}, len(ax)={}, len(rot)={}, len(circle)={}".format(
len(atom_permut),len(ax),len(rot),len(circle)))
#logger.info("len(circle)={}".format(len(circle))
#for aa, a, r, C in zip(atom_permut,ax,rot,circle):
for i, a in enumerate(ax):
# rotation strings see https://gitlab.com/ase/ase/blob/master/ase/utils/__init__.py#L235-261
plot_atoms(atom_permut[i],a,rotation=rot[i],radii=0.5,show_unit_cell=0,offset=(0,0))
for j, c in enumerate(circle):
logger.info("len(circle[{}])={}".format(j,len(c)))
a.add_patch(c[i])
for a,l,bb in zip(ax,label,plot_bb):
a.set_xlabel(l[0])
a.set_ylabel(l[1])
a.set_xlim(*bb[0,:])
a.set_ylim(*bb[1,:])
return fig, ax
def pack_sphere(C,
R_inner_constraint, # shell inner radius
R_outer_constraint, # shell outer radius
sfN, # number of surfactant molecules
inner_atom_number, # inner atom
outer_atom_number, # outer atom
surfactant = 'SDS',
counterion = 'NA',
tolerance = 2):
"""Creates context for filling Jinja2 PACKMOL input template in order to
generate preassembled surfactant spheres with couinterions at polar heads"""
logger = logging.getLogger(__name__)
logger.info(
"sphere with {:d} surfactant molecules in total.".format(sfN ) )
# sbX, sbY, sbZ = sb_measures
# spheres parallelt to x-axis
sphere = {}
ionsphere = {}
# surfactant spheres
# inner constraint radius: R + 1*tolerance
# outer constraint radius: R + 1*tolerance + l_surfactant
# ions between cylindric planes at
# inner radius: R + 1*tolerance + l_surfactant
# outer radius: R + 2*tolerance + l_surfactant
sphere["surfactant"] = surfactant
sphere["inner_atom_number"] = inner_atom_number
sphere["outer_atom_number"] = outer_atom_number
sphere["N"] = sfN
sphere["c"] = C
sphere["r_inner"] = R_inner
sphere["r_inner_constraint"] = R_inner_constraint
sphere["r_outer_constraint"] = R_outer_constraint
sphere["r_outer"] = R_outer
logging.info(
"sphere with {:d} molecules at {}, radius {}".format(
sphere["N"], sphere["c"], sphere["r_outer"]))
# ions at outer surface
ionsphere["ion"] = counterion
ionsphere["N"] = sphere["N"]
ionsphere["c"] = sphere["c"]
ionsphere["r_inner"] = sphere["r_outer"]
ionsphere["r_outer"] = sphere["r_outer"] + tolerance
# experience shows: movebadrandom advantegous for (hemi-) spheres
context = {
'spheres': [sphere],
'ionspheres': [ionsphere],
'movebadrandom': True,
}
return context
def memuse():
"""Quick overview on memory usage of objects in Jupyter notebook"""
# https://stackoverflow.com/questions/40993626/list-memory-usage-in-ipython-and-jupyter
# These are the usual ipython objects, including this one you are creating
ipython_vars = ['In', 'Out', 'exit', 'quit', 'get_ipython', 'ipython_vars']
# Get a sorted list of the objects and their sizes
return sorted([(x, sys.getsizeof(globals().get(x))) for x in dir(sys.modules['__main__']) if not x.startswith('_') and x not in sys.modules and x not in ipython_vars], key=lambda x: x[1], reverse=True)
os.environ['GMXLIB']
pmd.gromacs.GROMACS_TOPDIR = os.environ['GMXLIB']
prefix = '/mnt/dat/work/testuser/indenter/sandbox/20191110_packmol'
os.chdir(prefix)
hpc_max_specs = {
'forhlr2': {
'fw_queue_category': 'forhlr2_queue',
'fw_noqueue_category': 'forhlr2_noqueue',
'queue':'develop',
'physical_cores_per_node': 20,
'logical_cores_per_node': 40,
'nodes': 4,
'walltime': '00:60:00'
},
'juwels_devel': {
'fw_queue_category': 'juwels_queue',
'fw_noqueue_category': 'juwels_noqueue',
'queue':'devel',
'physical_cores_per_node': 48,
'logical_cores_per_node': 96,
'nodes': 8,
'walltime': '00:30:00'
},
'juwels': {
'fw_queue_category': 'juwels_queue',
'fw_noqueue_category': 'juwels_noqueue',
'queue':'batch',
'physical_cores_per_node': 48,
'logical_cores_per_node': 96,
'nodes': 1024,
'walltime': '00:30:00'
}
}
std_exports = {
'forhlr2': {
'OMP_NUM_THREADS': 1,
'KMP_AFFINITY': "'verbose,compact,1,0'",
'I_MPI_PIN_DOMAIN':'core'
},
'juwels': {
'OMP_NUM_THREADS': 1,
'KMP_AFFINITY': "'verbose,compact,1,0'",
'I_MPI_PIN_DOMAIN':'core'
}
}
# the FireWorks LaunchPad
lp = LaunchPad.auto_load() #Define the server and database
# FilePad behaves analogous to LaunchPad
fp = FilePad.auto_load()
The vollowing bash / tcl snippet converts a LAMMPS data file to PDB, assigning the desired names as mapped in a yaml file
#!/bin/bash
# echo "package require jlhvmd; jlh lmp2pdb indenter.lammps indenter.pdb" | vmd -eofexit
vmd -eofexit << 'EOF'
package require jlhvmd
topo readlammpsdata indenter.lammps
jlh type2name SDS_type2name.yaml
jlh name2res SDS_name2res.yaml
set sel [atomselect top all]
$sel writepdb indenter.pdb
EOF
pdb_chain.py indenter.pdb > indenter_wo_chainid.pdb
pdb_reres_by_atom_9999.py indenter_wo_chainid.pdb > indenter_reres.pdb
Requires
infile = os.path.join(prefix,'indenter_reres.pdb')
atoms = ase.io.read(infile,format='proteindatabank')
atoms
v = view(atoms,viewer='ngl')
v.view._remote_call("setSize", target="Widget", args=["400px", "400px"])
v.view.center()
v.view.background='#ffc'
v
S = atoms.get_positions()
C, R_sq = miniball.get_bounding_ball(S)
C # sphere center
R = np.sqrt(R_sq)
R # sphere radius
xmin = atoms.get_positions().min(axis=0)
xmax = atoms.get_positions().max(axis=0)
xmin
del S
# plot side views with sphere projections
plot_side_views_with_spheres(atoms,C,R)
# bounding sphere surface
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
us = np.array([
np.outer(np.cos(u), np.sin(v)),
np.outer(np.sin(u), np.sin(v)),
np.outer(np.ones(np.size(u)), np.cos(v))])
bs = C + R*us.T
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
x = atoms.get_positions()
ax.scatter(x[:,0], x[:,1],x[:,2], c='y', marker='o')
ax.plot_surface(*bs.T, color='b',alpha=0.1) # bounding sphere
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
tol = 2 # Ang
Utilize parmed to read pbd files ASE has difficulties to decipher.
infile = os.path.join(prefix,'1_SDS.pdb')
surfactant_pmd = pmd.load_file(infile)
surfactant_pmd.atoms[-1].atomic_number
surfactant_ase = ase.Atoms(
numbers=[1 if a.atomic_number == 0 else a.atomic_number for a in surfactant_pmd.atoms],
positions=surfactant_pmd.get_coordinates(0))
C_surfactant, R_sq_surfactant = miniball.get_bounding_ball(surfactant_ase.get_positions())
C_surfactant
R_surfactant = np.sqrt(R_sq_surfactant)
R_surfactant
C_surfactant
surfactant_ase[:5][1]
R_OSL = np.linalg.norm(C_surfactant - surfactant_ase[1].position)
R_OSL
d_head = R_surfactant - R_OSL # roughly: diameter of head group
R_inner = R + tol # place surfactant molecules outside of this sphere
R_inner_constraint = R + tol + d_head # place surfactant tail hydrocarbon within this sphere
R_outer_constraint = R + 2*R_surfactant + tol # place head group sulfur outside this sphere
R_outer = R + 2*R_surfactant + 2*tol # place suractant molecules within this sphere
rr = [R,R_inner,R_inner_constraint,R_outer_constraint,R_outer]
cc = [C]*5
plot_side_views_with_spheres(atoms,cc,rr,figsize=(20,8))
plt.show()
The template looks like this:
with open(os.path.join(infile_prefix,'surfactants_on_sphere.inp'),'r') as f:
print(f.read())
# get all placholders in template
template_file = os.path.join(infile_prefix,'surfactants_on_sphere.inp')
v = find_undeclared_variables(template_file)
v # we want to fill in these placeholder variables
surfactant = 'SDS'
counterion = 'NA'
tolerance = 2 # Ang
sfN = 200
l_surfactant = 2*R_surfactant
# head atom to be geometrically constrained
surfactant_head_bool_ndx = np.array([ a.name == 'S' for a in surfactant_pmd.atoms ],dtype=bool)
# tail atom to be geometrically constrained
surfactant_tail_bool_ndx = np.array([ a.name == 'C12' for a in surfactant_pmd.atoms ],dtype=bool)
head_atom_number = surfactant_head_ndx = np.argwhere(surfactant_head_bool_ndx)[0,0]
tail_atom_number = surfactant_tail_ndx = np.argwhere(surfactant_tail_bool_ndx)[0,0]
# settings can be overridden
packmol_script_context = {
'header': '20191113 TEST PACKING',
'system_name': '200_SDS_on_50_Ang_AFM_tip_model',
'tolerance': tolerance,
'write_restart': True,
'static_components': [
{
'name': 'indenter_reres'
}
]
}
# use pack_sphere function at the notebook's head to generate template context
packmol_script_context.update(
pack_sphere(
C,R_inner_constraint,R_outer_constraint, sfN,
tail_atom_number+1, head_atom_number+1, surfactant, counterion, tolerance))
packmol_script_context # context generated from system and constraint settings
env = jinja2.Environment()
template = jinja2.Template(open(template_file).read())
rendered = template.render(**packmol_script_context)
rendered_file = os.path.join(prefix,'rendered.inp')
with open(rendered_file,'w') as f:
f.write(rendered)
That's the rendered packmol input file:
print(rendered)
packmol = subprocess.Popen(['packmol'],
stdin=subprocess.PIPE,stdout=subprocess.PIPE, stderr=subprocess.PIPE,
cwd=prefix, encoding='utf-8')
outs, errs = packmol.communicate(input=rendered)
print(errs) # error with input from PIPE
packmol = subprocess.Popen(['packmol'],
stdin=open(rendered_file,'r'),stdout=subprocess.PIPE, stderr=subprocess.PIPE,
cwd=prefix, encoding='utf-8')
outs, errs = packmol.communicate(input=rendered)
print(outs)
with open('packmol.log','w') as f:
f.write(outs)
packmol_pdb = '200_SDS_on_50_Ang_AFM_tip_model_packmol.pdb'
infile = os.path.join(prefix, packmol_pdb)
surfactant_shell_pmd = pmd.load_file(infile)
# with ParmEd and nglview we get automatic bond guessing
pmd_view = nglview.show_parmed(surfactant_shell_pmd)
pmd_view.clear_representations()
pmd_view.background = 'white'
pmd_view.add_representation('ball+stick')
pmd_view
surfactant_shell_ase = ase.Atoms(
numbers=[1 if a.atomic_number == 0 else a.atomic_number for a in surfactant_shell_pmd.atoms],
positions=surfactant_shell_pmd.get_coordinates(0))
# with ASE, we get no bonds at all
ase_view = nglview.show_ase(surfactant_shell_ase)
ase_view.clear_representations()
ase_view.background = 'white'
ase_view.add_representation('ball+stick')
ase_view
Get bounding sphere again and display AFM tip bounding spphere as well as surfactant layer bounding sphere
C_shell, R_sq_shell = miniball.get_bounding_ball(surfactant_shell_ase.get_positions())
C_shell
R_shell = np.sqrt(R_sq_shell)
R_shell
plot_side_views_with_spheres(surfactant_shell_ase,[C,C_shell],[R,R_shell])
surfactant_shell_pmd
R # Angstrom
A_Ang = 4*np.pi*R**2 # area in Ansgtrom
A_nm = A_Ang / 10**2
A_nm
n_per_nm_sq = np.array([0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0]) # molecules per square nm
N = np.round(A_nm*n_per_nm_sq).astype(int)
N # molecule numbers corresponding to surface concentrations
project_id = 'juwels-packmol-2020-03-09'
# queries to the data base are simple dictionaries
query = {
'metadata.project': project_id,
}
# use underlying MongoDB functionality to check total number of documents matching query
fp.filepad.count_documents(query)
infile_prefix = os.path.join(prefix,'packmol_infiles')
infiles = sorted(glob.glob(os.path.join(infile_prefix,'*.inp')))
files = { os.path.basename(f): f for f in infiles }
# metadata common to all these files
metadata = {
'project': project_id,
'type': 'template'
}
fp_files = []
# insert these input files into data base
for name, file_path in files.items():
identifier = '/'.join((project_id,name)) # identifier is like a path on a file system
metadata["name"] = name
fp_files.append( fp.add_file(file_path,identifier=identifier,metadata = metadata) )
# queries to the data base are simple dictionaries
query = {
'metadata.project': project_id,
'metadata.type': 'template'
}
# use underlying MongoDB functionality to check total number of documents matching query
fp.filepad.count_documents(query)
print(identifier)
# on a lower level, each object has a unique "GridFS id":
pprint(fp_files) # underlying GridFS id and readable identifiers
# settings can be overridden
for n in N:
packmol_script_context = {
'header': '{:s} packing SDS around AFM probe model'.format(project_id),
'system_name': '{:d}_SDS_on_50_Ang_AFM_tip_model'.format(n),
'tolerance': tolerance,
'write_restart': True,
'static_components': [
{
'name': 'indenter_reres'
}
]
}
# use pack_sphere function at the notebook's head to generate template context
packmol_script_context.update(
pack_sphere(
C,R_inner_constraint,R_outer_constraint, n,
tail_atom_number, head_atom_number, surfactant, counterion, tolerance))
data_prefix = os.path.join(prefix,'packmol_datafiles')
datafiles = sorted(glob.glob(os.path.join(data_prefix,'*')))
files = { os.path.basename(f): f for f in datafiles }
# metadata common to all these files
metadata = {
'project': project_id,
'type': 'data'
}
fp_files = []
# insert these input files into data base
for name, file_path in files.items():
identifier = '/'.join((project_id,name)) # identifier is like a path on a file system
metadata["name"] = name
fp_files.append( fp.add_file(file_path,identifier=identifier,metadata = metadata) )
fp_files
machine = 'juwels_devel'
parametric_dimension_labels = ['nmolecules']
parametric_dimensions = [ {
'nmolecules': N } ]
# for testing
parametric_dimensions = [ {
'nmolecules': [N[0]] } ]
parameter_sets = list(
itertools.chain(*[
itertools.product(*list(
p.values())) for p in parametric_dimensions ]) )
parameter_dict_sets = [ dict(zip(parametric_dimension_labels,s)) for s in parameter_sets ]
wf_name = 'PACKMOL {machine:}, {id:}'.format(machine=machine,id=project_id)
fw_name_template = 'nmolecules: {nmolecules:d}'
fw_list = []
fts = [
GetFilesByQueryTask(
query = {
'metadata->project': project_id,
'metadata->name': 'surfactants_on_sphere.inp'
},
limit = 1,
new_file_names = ['input.template'] ),
GetFilesByQueryTask(
query = {
'metadata->project': project_id,
'metadata->name': 'indenter_reres.pdb'
},
limit = 1,
new_file_names = ['indenter.pdb'] ),
GetFilesByQueryTask(
query = {
'metadata->project': project_id,
'metadata->name': '1_SDS.pdb'
},
limit = 1,
new_file_names = ['1_SDS.pdb'] ),
GetFilesByQueryTask(
query = {
'metadata->project': project_id,
'metadata->name': '1_NA.pdb'
},
limit = 1,
new_file_names = ['1_NA.pdb'] )
]
files_out = {
'input_file': 'input.template',
'indenter_file': 'indenter.pdb',
'surfatcant_file': '1_SDS.pdb',
'counterion_file': '1_NA.pdb'}
fw_root = Firework(fts,
name = wf_name,
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'input_file_query'
}
}
)
fw_list.append(fw_root)
## Parametric sweep
for d in parameter_dict_sets:
### Template
files_in = {'input_file': 'input.template' }
files_out = { 'input_file': 'input.inp' }
# exports = std_exports[machine].copy()
# Jinja2 context:
packmol_script_context = {
'header': '{:s} packing SDS around AFM probe model'.format(project_id),
'system_name': '{:d}_SDS_on_50_Ang_AFM_tip_model'.format(n),
'tolerance': tolerance,
'write_restart': True,
'static_components': [
{
'name': 'indenter'
}
]
}
# use pack_sphere function at the notebook's head to generate template context
packmol_script_context.update(
pack_sphere(
C,R_inner_constraint,R_outer_constraint, d["nmolecules"],
tail_atom_number+1, head_atom_number+1, surfactant, counterion, tolerance))
ft_template = TemplateWriterTask( {
'context': packmol_script_context,
'template_file': 'input.template',
'template_dir': '.',
'output_file': 'input.inp'} )
fw_template = Firework([ft_template],
name = ', '.join(('template', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'fill_template',
**d
}
},
parents = [ fw_root ] )
fw_list.append(fw_template)
### PACKMOL
files_in = {
'input_file': 'input.inp',
'indenter_file': 'indenter.pdb',
'surfatcant_file': '1_SDS.pdb',
'counterion_file': '1_NA.pdb' }
files_out = {
'data_file': '*_packmol.pdb'}
ft_pack = CmdTask(
cmd='packmol',
opt=['< input.inp'],
stderr_file = 'std.err',
stdout_file = 'std.out',
store_stdout = True,
store_stderr = True,
use_shell = True,
fizzle_bad_rc= True)
fw_pack = Firework([ft_pack],
name = ', '.join(('packmol', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_queue_category'],
'_queueadapter': {
'queue': hpc_max_specs[machine]['queue'],
'walltime' : hpc_max_specs[machine]['walltime'],
'ntasks': 1,
},
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'packmol',
**d
}
},
parents = [ fw_root, fw_template] )
fw_list.append(fw_pack)
### Store
files_in = {'data_file': 'packed.pdb' }
ft_transfer = AddFilesTask( {
'compress': True ,
'paths': "packed.pdb",
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'type': 'initial_config',
**d }
} )
fw_transfer = Firework([ft_transfer],
name = ', '.join(('transfer', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'transfer',
**d
}
},
parents = [ fw_pack ] )
fw_list.append(fw_transfer)
wf = Workflow(fw_list,
name = wf_name,
metadata = {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'type': 'packing'
})
wf.to_file('packing.json')
lp.add_wf(wf)
query = {
"metadata.project": project_id,
}
fp.filepad.count_documents(query)
query = {
"metadata.project": project_id,
"metadata.type": 'initial_config',
}
fp.filepad.count_documents(query)
parameter_names = ['nmolecules']
surfactant_shell_pmd_list = []
match_aggregation = {
"$match": query
}
sort_aggregation = {
"$sort": {
"metadata.datetime": pymongo.DESCENDING,
}
}
group_aggregation = {
"$group": {
"_id": { p: '$metadata.{}'.format(p) for p in parameter_names },
"degeneracy": {"$sum": 1}, # number matching data sets
"latest": {"$first": "$gfs_id"} # unique gridfs id of file
}
}
aggregation_pipeline = [ match_aggregation, sort_aggregation, group_aggregation ]
cursor = fp.filepad.aggregate(aggregation_pipeline)
for i, c in enumerate(cursor):
content, metadata = fp.get_file_by_id(c["latest"])
with tempfile.NamedTemporaryFile() as tmp:
tmp.write(content)
surfactant_shell_pmd_list.append(pmd.load_file(tmp.name))
print('.',end='')
print('')
system_selection = surfactant_shell_pmd_list
nx = 3
ny = len(system_selection)
view_labels = ['xy','xz','yz']
molecule_counter = lambda s: np.count_nonzero([r.name == 'SDS' for r in s.residues])
label_generator = lambda i,j: '{:d} SDS, {:s} projection'.format(
molecule_counter(system_selection[i]),view_labels[j])
figsize = (4*nx,4*ny)
fig, axes = plt.subplots(ny,3,figsize=figsize)
for i,system in enumerate(system_selection):
system_ase = ase.Atoms(
numbers=[1 if a.atomic_number == 0 else a.atomic_number for a in system.atoms],
positions=system.get_coordinates(0))
C_shell, R_sq_shell = miniball.get_bounding_ball(system_ase.get_positions())
R_shell = np.sqrt(R_sq_shell)
plot_side_views_with_spheres(
system_ase,[C,C_shell],[R,R_shell],fig=fig,ax=axes[i,:])
for j, ax in enumerate(axes[i,:]):
ax.set_title(label_generator(i,j))
del system_ase
gc.collect()
del fig
del axes
gromacs.config.logfilename
gromacs.environment.flags
# if true, then stdout and stderr are returned as strings by gromacs wrapper commands
gromacs.environment.flags['capture_output'] = False
print(gromacs.release())
prefix
system = '200_SDS_on_50_Ang_AFM_tip_model'
pdb = system + '.pdb'
gro = system + '.gro'
top = system + '.top'
posre = system + '.posre.itp'
# Remove any chain ID from pdb and tidy up
pdb_chain = subprocess.Popen(['pdb_chain',],
stdin=open(packmol_pdb,'r'),stdout=subprocess.PIPE, stderr=subprocess.PIPE,
cwd=prefix, encoding='utf-8')
pdb_tidy = subprocess.Popen(['pdb_tidy',],
stdin=pdb_chain.stdout,stdout=open(pdb,'w'), stderr=subprocess.PIPE,
cwd=prefix, encoding='utf-8')
rc,out,err=gromacs.pdb2gmx(
f=pdb,o=gro,p=top,i=posre,ff='charmm36',water='tip3p',
stdout=False,stderr=False)
print(out)
gro_boxed = system + '_boxed.gro'
rc,out,err=gromacs.editconf(
f=gro,o=gro_boxed,d=2.0,bt='cubic',
stdout=False,stderr=False)
print(out)
machine = 'juwels_devel'
parametric_dimension_labels = ['nmolecules']
parametric_dimensions = [ {
'nmolecules': N } ]
# for testing
parametric_dimensions = [ {
'nmolecules': [N[0]] } ]
parameter_sets = list(
itertools.chain(*[
itertools.product(*list(
p.values())) for p in parametric_dimensions ]) )
parameter_dict_sets = [ dict(zip(parametric_dimension_labels,s)) for s in parameter_sets ]
source_project_id = 'juwels-packmol-2020-03-09'
project_id = 'juwels-gromacs-prep-2020-03-11'
wf_name = 'GROMACS preparations {machine:}, {id:}'.format(machine=machine,id=project_id)
fw_name_template = 'nmolecules: {nmolecules:d}'
fw_list = []
fts = [
CmdTask(
cmd='echo',
opt=['"Dummy root"'],
store_stdout = False,
store_stderr = False,
use_shell = True,
fizzle_bad_rc= True) ]
files_out = []
fw_root = Firework(fts,
name = wf_name,
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'dummy_root'
}
}
)
fw_list.append(fw_root)
## Parametric sweep
for d in parameter_dict_sets:
### File retrieval
#files_in = {'input_file': 'input.template' }
files_in = {}
files_out = { 'data_file': 'in.pdb' }
# exports = std_exports[machine].copy()
fts_fetch = [ GetFilesByQueryTask(
query = {
'metadata->project': source_project_id,
'metadata->type': 'initial_config',
'metadata->nmolecules': d["nmolecules"]
},
sort_key = 'metadata.datetime',
sort_direction = pymongo.DESCENDING,
limit = 1,
new_file_names = ['in.pdb'] ) ]
fw_fetch = Firework(fts_fetch,
name = ', '.join(('fetch', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'fetch',
**d
}
},
parents = [ fw_root ] )
fw_list.append(fw_fetch)
### PDB chain
files_in = {'data_file': 'in.pdb' }
files_out = {'data_file': 'out.pdb'}
fts_pdb_chain = CmdTask(
cmd='pdb_chain',
opt=['< in.pdb > out.pdb'],
store_stdout = False,
store_stderr = False,
use_shell = True,
fizzle_bad_rc= True)
fw_pdb_chain = Firework(fts_pdb_chain,
name = ', '.join(('pdb_chain', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'pdb_chain',
**d
}
},
parents = [ fw_fetch ] )
fw_list.append(fw_pdb_chain)
### PDB tidy
files_in = {'data_file': 'in.pdb' }
files_out = {'data_file': 'out.pdb'}
fts_pdb_tidy = CmdTask(
cmd='pdb_tidy',
opt=['< in.pdb > out.pdb'],
store_stdout = False,
store_stderr = False,
use_shell = True,
fizzle_bad_rc= True)
fw_pdb_tidy = Firework(fts_pdb_tidy,
name = ', '.join(('pdb_tidy', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'pdb_tidy',
**d
}
},
parents = [ fw_pdb_chain ] )
fw_list.append(fw_pdb_tidy)
### GMX pdb2gro
files_in = {'data_file': 'in.pdb' }
files_out = {
'coordinate_file': 'out.gro',
'topology_file': 'out.top',
'restraint_file': 'out.posre.itp'}
fts_gmx_pdb2gro = [ CmdTask(
cmd='gmx',
opt=['pdb2gmx',
'-f', 'in.pdb',
'-o', 'out.gro',
'-p', 'out.top',
'-i', 'out.posre.itp',
'-ff', 'charmm36',
'-water' , 'tip3p'],
stderr_file = 'std.err',
stdout_file = 'std.out',
store_stdout = True,
store_stderr = True,
use_shell = True,
fizzle_bad_rc= True) ]
fw_gmx_pdb2gro = Firework(fts_gmx_pdb2gro,
name = ', '.join(('gmx_pdb2gro', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'gmx_pdb2gro',
**d
}
},
parents = [ fw_pdb_tidy ] )
fw_list.append(fw_gmx_pdb2gro)
### GMX editconf
files_in = {
'coordinate_file': 'in.gro',
'topology_file': 'in.top',
'restraint_file': 'in.posre.itp'}
files_out = {
'coordinate_file': 'out.gro',
'topology_file': 'in.top',
'restraint_file': 'in.posre.itp'}
fts_gmx_editconf = [ CmdTask(
cmd='gmx',
opt=['editconf',
'-f', 'in.gro',
'-o', 'out.gro',
'-d', 2.0, # distance between content and box boundary in nm
'-bt', 'cubic', # box type
],
stderr_file = 'std.err',
stdout_file = 'std.out',
store_stdout = True,
store_stderr = True,
use_shell = True,
fizzle_bad_rc= True) ]
fw_gmx_editconf = Firework(fts_gmx_editconf,
name = ', '.join(('gmx_editconf', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'_files_out': files_out,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'gmx_editconf',
**d
}
},
parents = [ fw_gmx_pdb2gro ] )
fw_list.append(fw_gmx_editconf)
### Store
#files_in = {'data_file': 'packed.pdb' }
files_in = {
'coordinate_file': 'default.gro',
'topology_file': 'default.top',
'restraint_file': 'default.posre.itp' }
fts_push = [
AddFilesTask( {
'compress': True ,
'paths': "default.gro",
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'type': 'initial_config_gro',
**d }
} ),
AddFilesTask( {
'compress': True ,
'paths': "default.top",
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'type': 'initial_config_top',
**d }
} ),
AddFilesTask( {
'compress': True ,
'paths': "default.posre.itp",
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'type': 'initial_config_posre_itp',
**d }
} ) ]
fw_push = Firework(fts_push,
name = ', '.join(('push', fw_name_template.format(**d))),
spec = {
'_category': hpc_max_specs[machine]['fw_noqueue_category'],
'_files_in': files_in,
'metadata': {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'step': 'push',
**d
}
},
parents = [ fw_gmx_editconf ] )
fw_list.append(fw_push)
wf = Workflow(fw_list,
name = wf_name,
metadata = {
'project': project_id,
'datetime': str(datetime.datetime.now()),
'type': 'gmx_prep'
})
wf.as_dict()
lp.add_wf(wf)
wf.to_file('wf-{:s}.yaml'.format(project_id))
query = {
"metadata.project": project_id,
}
fp.filepad.count_documents(query)
query = {
"metadata.project": project_id,
"metadata.type": 'initial_config_gro',
}
fp.filepad.count_documents(query)
parameter_names = ['nmolecules']
gro_list = []
match_aggregation = {
"$match": query
}
sort_aggregation = {
"$sort": {
"metadata.datetime": pymongo.DESCENDING,
}
}
group_aggregation = {
"$group": {
"_id": { p: '$metadata.{}'.format(p) for p in parameter_names },
"degeneracy": {"$sum": 1}, # number matching data sets
"latest": {"$first": "$gfs_id"} # unique gridfs id of file
}
}
aggregation_pipeline = [ match_aggregation, sort_aggregation, group_aggregation ]
cursor = fp.filepad.aggregate(aggregation_pipeline)
for i, c in enumerate(cursor):
content, metadata = fp.get_file_by_id(c["latest"])
with tempfile.NamedTemporaryFile(suffix='.gro') as tmp:
tmp.write(content)
gro_list.append(pmd.load_file(tmp.name))
print('.',end='')
print('')
gro_list
# with ParmEd and nglview we get automatic bond guessing
pmd_view = nglview.show_parmed(gro_list[0])
pmd_view.clear_representations()
pmd_view.background = 'white'
pmd_view.add_representation('ball+stick')
pmd_view
system_selection = gro_list
nx = 3
ny = len(system_selection)
view_labels = ['xy','xz','yz']
molecule_counter = lambda s: np.count_nonzero([r.name == 'SDS' for r in s.residues])
label_generator = lambda i,j: '{:d} SDS, {:s} projection'.format(
molecule_counter(system_selection[i]),view_labels[j])
figsize = (4*nx,4*ny)
fig, axes = plt.subplots(ny,3,figsize=figsize)
for i,system in enumerate(system_selection):
system_ase = ase.Atoms(
numbers=[1 if a.atomic_number == 0 else a.atomic_number for a in system.atoms],
positions=system.get_coordinates(0))
C_shell, R_sq_shell = miniball.get_bounding_ball(system_ase.get_positions())
R_shell = np.sqrt(R_sq_shell)
plot_side_views_with_spheres(
system_ase,[C,C_shell],[R,R_shell],fig=fig,ax=axes[i,:])
for j, ax in enumerate(axes[i,:]):
ax.set_title(label_generator(i,j))
del system_ase
gc.collect()
Just to be safe, relax the system a little with positional constraints applied to all ions.
os.getcwd()
em_mdp = gromacs.fileformats.MDP('em.mdp.template')
# no change
em_mdp.write('em.mdp')
gmx_grompp = gromacs.grompp.Popen(
f='em.mdp',c=gro,r=gro,o='em.tpr',p=top,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_grompp.stdout.read()
err = gmx_grompp.stderr.read()
print(err)
print(out)
gmx_mdrun = gromacs.mdrun.Popen(
deffnm='em',v=True,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
for line in gmx_mdrun.stdout:
print(line.decode(), end='')
out = gmx_mdrun.stdout.read()
err = gmx_mdrun.stderr.read()
print(err)
print(out)
em_file = 'em.edr'
em_df = panedr.edr_to_df(em_file)
em_df.columns
fig, ax = plt.subplots(3,2,figsize=(10,12))
em_df.plot('Time','Potential',ax=ax[0,0])
em_df.plot('Time','Pressure',ax=ax[0,1])
em_df.plot('Time','Bond',ax=ax[1,0])
em_df.plot('Time','Position Rest.',ax=ax[1,1])
#em_df.plot('Time','COM Pull En.',ax=ax[1,1])
em_df.plot('Time','Coulomb (SR)',ax=ax[2,0])
em_df.plot('Time','Coul. recip.',ax=ax[2,1])
mda_trr = mda.Universe(gro,'em.trr')
mda_view = nglview.show_mdanalysis(mda_trr)
mda_view.clear_representations()
mda_view.background = 'white'
mda_view.add_representation('ball+stick')
mda_view
Utilize harmonic pulling to attach surfactants to substrate closely.
#pdb = '200_SDS_on_50_Ang_AFM_tip_model.pdb'
gro = 'em.gro'
top = 'sys.top'
ndx = 'standard.ndx'
pmd_top_gro = pmd.gromacs.GromacsTopologyFile(top)
#pmd_top_pdb = pmd.gromacs.GromacsTopologyFile(top)
pmd_gro = pmd.gromacs.GromacsGroFile.parse(gro)
pmd_top_gro.box = pmd_gro.box
pmd_top_gro.positions = pmd_gro.positions
#pmd_pdb = pmd.formats.pdb.PDBFile.parse(pdb)
#pmd_top_pdb.box = pmd_pdb.box
#pmd_top_pdb.positions = pmd_pdb.positions
tail_atom_ndx = np.array([
i+1 for i,a in enumerate(pmd_top_gro.atoms) if a.name == 'C12' and a.residue.name == 'SDS' ])
# gromacs ndx starts at 1
tail_atom_ndx
gmx_make_ndx = gromacs.make_ndx.Popen(
f=gro,o=ndx,
input='q',
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out_str, err_str = gmx_make_ndx.communicate()
print(out_str)
print(err_str)
pull_groups_ndx_in = gromacs.fileformats.NDX(ndx)
pull_groups_ndx_out = gromacs.fileformats.NDX()
for i, a in enumerate(tail_atom_ndx):
pull_group_name = 'pull_group_{:04d}'.format(i)
pull_groups_ndx_out[pull_group_name] = a
pull_groups_ndx_in.update(pull_groups_ndx_out)
pull_groups_ndx_in.write('pull_groups.ndx')
# gromacs wrapper parses mdp files
pull_mdp = gromacs.fileformats.MDP('pull.mdp.template')
pull_mdp['nsteps'] = 10000
N_pull_coords = len(pull_groups_ndx_out)
pull_mdp['pull_ncoords'] = N_pull_coords
pull_mdp['pull_ngroups'] = N_pull_coords + 1
pull_mdp['pull_group1_name'] = 'Substrate' # the reference group
for i, n in enumerate(pull_groups_ndx_out):
pull_mdp["pull_group{:d}_name".format(i+2)] = n
pull_mdp["pull_coord{:d}_type".format(i+1)] = 'umbrella' # harmonic potential
pull_mdp["pull_coord{:d}_geometry".format(i+1)] = 'distance' # simple distance increase
pull_mdp["pull_coord{:d}_dim".format(i+1)] = 'Y Y Y' # pull in all directions
pull_mdp["pull_coord{:d}_groups".format(i+1)] = "1 {:d}".format(i+2) # groups 1 (Chain A) and 2 (Chain B) define the reaction coordinate
pull_mdp["pull_coord{:d}_start".format(i+1)] = 'yes' # define initial COM distance > 0
pull_mdp["pull_coord{:d}_rate".format(i+1)] = -0.1 # 0.1 nm per ps = 10 nm per ns
pull_mdp["pull_coord{:d}_k".format(i+1)] = 1000 # kJ mol^-1 nm^-2
pull_mdp.write('pull.mdp')
gmx_grompp = gromacs.grompp.Popen(
f='pull.mdp',n='pull_groups.ndx',c=gro,r=gro,o='pull.tpr',p=top,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_grompp.stdout.read()
err = gmx_grompp.stderr.read()
print(err)
print(out)
gmx_mdrun = gromacs.mdrun.Popen(
deffnm='pull',v=True,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_mdrun.stdout.read()
err = gmx_mdrun.stderr.read()
print(err)
print(out)
edr_file = 'pull.edr'
edr_df = panedr.edr_to_df(edr_file)
edr_df.columns
fig, ax = plt.subplots(3,2,figsize=(10,12))
edr_df.plot('Time','Potential',ax=ax[0,0])
edr_df.plot('Time','Pressure',ax=ax[0,1])
#edr_df.plot('Time','Bond',ax=ax[1,0])
edr_df.plot('Time','Position Rest.',ax=ax[1,0])
edr_df.plot('Time','COM Pull En.',ax=ax[1,1])
edr_df.plot('Time','Coulomb (SR)',ax=ax[2,0])
edr_df.plot('Time','Coul. recip.',ax=ax[2,1])
# read xvg file
pull_f_xvg = mda.auxiliary.XVG.XVGFileReader('pull_pullf.xvg')
pull_f_t = pull_f_xvg.read_all_times()
# first data column contains time, strip
pull_f = np.array([ f.data[1:] for f in pull_f_xvg ])
for i in range(0,199,50):
plt.plot(pull_f_t,pull_f[:,i])
pull_x_xvg = gromacs.fileformats.XVG('pull_pullx.xvg',)
pull_x_xvg.array
len(pull_x_xvg.names)
# that many columns perr pull coordinate
N_cols_per_coord = int(len(pull_x_xvg.names) / N_pull_coords)
# with content
legend = pull_x_xvg.names[:11]
legend
pull_x_xvg.names[-3:]
for i in range(11):
plt.plot(pull_x_xvg.array[0,:],pull_x_xvg.array[i+1,:],label=legend[i])
plt.legend()
for i in range(11):
plt.plot(pull_x_xvg.array[0,:],pull_x_xvg.array[i+1,:],label=legend[i])
plt.legend()
# sqrt(dx^2+dy^2+dz^2), the distance between pulling groups (i.e. one surfactant tail atom and the Au COM)
for i in range(0,150,30):
plt.plot(pull_x_xvg.array[0,:],
np.sqrt(pull_x_xvg.array[i*11+3,:]**2+pull_x_xvg.array[i*11+4,:]**2+pull_x_xvg.array[i*11+5,:]**2))
plt.legend()
pull_x_xvg.plot(columns=[0,12])
gro_em = 'pull.gro'
mda_trr = mda.Universe('em.gro','pull.trr')
mda_view = nglview.show_mdanalysis(mda_trr)
mda_view.clear_representations()
mda_view.background = 'white'
mda_view.add_representation('ball+stick')
mda_view
mda_xtc = mda.Universe(gro,'pull.xtc')
mda_view = nglview.show_mdanalysis(mda_xtc)
mda_view.clear_representations()
mda_view.background = 'white'
mda_view.add_representation('ball+stick')
mda_view
substrate = mda_trr.atoms[mda_trr.atoms.names == 'AU']
surfactant_head = mda_trr.atoms[mda_trr.atoms.names == 'S']
rms_substrate = mda_rms.RMSD(substrate,ref_frame=0)
rms_substrate.run()
rmsd = rms_substrate.rmsd.T # transpose makes it easier for plotting
time = rmsd[1]
plt.plot(time,rmsd[2])
rms_surfactant_head = mda_rms.RMSD(surfactant_head,ref_frame=0)
rms_surfactant_head.run()
rmsd = rms_surfactant_head.rmsd.T # transpose makes it easier for plotting
time = rmsd[1]
plt.plot(time,rmsd[2])
len(mda_trr.trajectory)
rdf_substrate_headgroup = mda_rdf.InterRDF(
substrate,surfactant_head,range=(0.0,80.0),verbose=True)
bins = []
rdf = []
for i in range(len(mda_trr.trajectory)):
rdf_substrate_headgroup = mda_rdf.InterRDF(
substrate,surfactant_head,range=(0.0,80.0),verbose=True)
rdf_substrate_headgroup.run(start=i,stop=i+1)
bins.append(rdf_substrate_headgroup.bins.copy())
rdf.append(rdf_substrate_headgroup.rdf.copy())
bins = np.array(bins)
rdf = np.array(rdf)
# indicates desired approach towards substrate
plt.plot(bins[0],rdf[0],label="Initial RDF")
plt.plot(bins[3],rdf[4],label="Intermediat RDF")
plt.plot(bins[-1],rdf[-1],label="Final RDF")
plt.legend()
Now, fill the box with water.
gro = 'pull.gro'
# use -scale 0.5 -maxsol N for non-standard conditions
gmx_solvate = gromacs.solvate.Popen(
cp=gro, cs='spc216.gro',o='solvated.gro',p=top,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_solvate.stdout.read()
err = gmx_solvate.stderr.read()
print(out)
print(err)
Again, relax the system a little with positional constraints applied to all ions.
lpad = LaunchPad.auto_load()
A trial task sent to FORHLR2:
gmx_test_task = CmdTask(
cmd = 'gmx',
opt = '-h',
stderr_file = 'std.err',
stdout_file = 'std.out',
use_shell = True)
gmx_test_fw = Firework(
gmx_test_task,
name = 'FORHLR2 GMX test fw',
spec = {
'_category': 'forhlr2_noqueue',
'_files_out': {
'stdout': 'std.out',
'stderr': 'std.err'}
} )
fw_ids = lpad.add_wf(gmx_test_fw)
fw_ids
# lpad.delete_wf(INSERT_ID,delete_launch_dirs=True)
top = 'sys.top'
gro = 'solvated.gro'
gmx_grompp = gromacs.grompp.Popen(
f='em_solvated.mdp',c=gro,r=gro,o='em_solvated.tpr',p=top,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_grompp.stdout.read()
err = gmx_grompp.stderr.read()
print(err)
print(out)
Utilize fabric to transfer files files to remote resource conveniently:
c = fabric.Connection('forhlr2') # host defined in ssh config
res = c.run('ws_find fw') # get remote directory of Fireworks workspace
res.command
now = datetime.now().isoformat()
remote_path = os.path.sep.join((res.stdout.strip(),'file_transfer',now))
remote_path
res = c.run(' '.join(['mkdir','-p',remote_path]))
file_name = 'em_solvated.tpr'
local_file = os.path.sep.join((prefix,file_name))
remote_file = os.path.sep.join((remote_path,file_name))
res = c.put(local_file,remote_file)
res.local
# FileTransferTask does not work anymore
#ft = FileTransferTask(
# mode = 'rtransfer',
# files = [ {'src':local_file, 'dest':remote_path} ],
# server = c.host,
# user = c.user )
#fw = Firework(
# ft,
# name = 'BWCLOUD remote transef to FORHLR2',
# spec = {
# '_category': 'bwcloud_std',
# } )
fw_ids = lpad.add_wf(fw)
ft = FileTransferTask(
mode = 'copy',
files = [ {'src':remote_file, 'dest':os.path.curdir} ] )
gmx_mdrun_task = CmdTask(
cmd = 'gmx',
opt = ['mdrun','-v','-deffnm','em_solvated'],
stderr_file = 'std.err',
stdout_file = 'std.out',
use_shell = True)
gmx_log_tracker = Tracker('em_solvated.log')
gmx_mdrun_fw = Firework(
[ft,gmx_mdrun_task],
name = 'FORHLR2 GMX mdrun em_solvated',
spec = {
'_category': 'forhlr2_queue',
'_queueadapter': {
'cpus_per_task': 1,
'ntasks_per_node': 20,
'ntasks': 40,
'queue': 'normal',
'walltime': '24:00'
},
'_files_out': {
'log': '*.log',
'trr': '*.trr',
'edr': '*.edr',
'gro': '*.gro' },
'_trackers' : [ gmx_log_tracker ]
} )
pprint(gmx_mdrun_fw.as_dict())
fw_ids = lpad.add_wf(gmx_mdrun_fw)
fw_id = list(fw_ids.values())[0]
instead of relying on the returned fw_id, we can also query the Firework added latest
fw_ids = lpad.get_fw_ids(sort=[('created_on',pymongo.DESCENDING)],limit=1)
fw_id = fw_ids[0]
fw_id
We query the remote directory our FireWork ran in
launch_dir = lpad.get_launchdir(fw_id)
launch_dir
c = fabric.Connection('forhlr2') # host defined in ssh config
res = c.run('ls -lht {}'.format(launch_dir)) # look at remote directory contents
glob_pattern = os.path.join(launch_dir,'em_solvated.*')
res = c.run('ls {}'.format(glob_pattern))
res.stdout
for f in res.stdout.splitlines():
c.get(f)
em_file = 'em_solvated.edr'
em_df = panedr.edr_to_df(em_file)
em_df.columns
fig, ax = plt.subplots(3,2,figsize=(10,12))
em_df.plot('Time','Potential',ax=ax[0,0])
em_df.plot('Time','Pressure',ax=ax[0,1])
em_df.plot('Time','Bond',ax=ax[1,0])
em_df.plot('Time','Position Rest.',ax=ax[1,1])
#em_df.plot('Time','COM Pull En.',ax=ax[1,1])
em_df.plot('Time','Coulomb (SR)',ax=ax[2,0])
em_df.plot('Time','Coul. recip.',ax=ax[2,1])
try:
del em_df
except:
pass
mda_trr = mda.Universe('solvated.gro','em_solvated.trr')
# check unique resiude names in system
resnames = np.unique([ r.resname for r in mda_trr.residues ])
resnames
mda_view = nglview.show_mdanalysis(mda_trr)
mda_view.clear_representations()
mda_view.background = 'white'
mda_view.add_representation(repr_type='ball+stick',selection='SDS')
mda_view.add_representation(repr_type='ball+stick',selection='NA')
mda_view.add_representation(repr_type='spacefill',selection='AUM',color='yellow')
mda_view.center()
mda_view
try:
del mda_trr
except:
pass
try:
del mda_view
except:
pass
top = 'sys.top'
gro = 'em_solvated.gro'
ndx = 'nvt.ndx'
lpad = LaunchPad.auto_load()
pmd_top_gro = pmd.gromacs.GromacsTopologyFile(top)
pmd_gro = pmd.gromacs.GromacsGroFile.parse(gro)
pmd_top_gro.box = pmd_gro.box
pmd_top_gro.positions = pmd_gro.positions
non_substrate_ndx = np.array([
i+1 for i,a in enumerate(pmd_top_gro.atoms) if a.residue.name != 'AUM' ])
# gromacs ndx starts at 1
len(pmd_top_gro.atoms)
len(non_substrate_ndx)
len(pmd_top_gro.atoms) - len(non_substrate_ndx) # double-check non-substrate and substrate atom numbers
try:
del pmd_top_gro
except:
pass
try:
del pmd_gro
except:
pass
gmx_make_ndx = gromacs.make_ndx.Popen(
f=gro,o=ndx,
input='q',
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out_str, err_str = gmx_make_ndx.communicate()
print(out_str)
print(err_str)
ndx_in = gromacs.fileformats.NDX(ndx)
ndx_in['non-Substrate'] = non_substrate_ndx
ndx_in.write(ndx)
try:
del ndx_in
except:
pass
try:
del non_substrate_ndx
except:
pass
gmx_grompp = gromacs.grompp.Popen(
f='nvt.mdp',n=ndx,c=gro,r=gro,o='nvt.tpr',p=top,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_grompp.stdout.read()
err = gmx_grompp.stderr.read()
print(err)
print(out)
Utilize fabric to transfer files files to remote resource conveniently:
c = fabric.Connection('forhlr2') # host defined in ssh config
res = c.run('ws_find fw') # get remote directory of Fireworks workspace
res.command
now = datetime.now().isoformat()
remote_path = os.path.sep.join((res.stdout.strip(),'file_transfer',now))
remote_path
res = c.run(' '.join(['mkdir','-p',remote_path]))
file_name = 'nvt.tpr'
local_file = os.path.sep.join((prefix,file_name))
remote_file = os.path.sep.join((remote_path,file_name))
res = c.put(local_file,remote_file)
res.local
ft = FileTransferTask(
mode = 'copy',
files = [ {'src':remote_file, 'dest':os.path.curdir} ] )
gmx_mdrun_task = CmdTask(
cmd = 'gmx',
opt = ['mdrun','-v','-deffnm','nvt'],
stderr_file = 'std.err',
stdout_file = 'std.out',
use_shell = True)
gmx_log_tracker = Tracker('nvt.log')
gmx_mdrun_fw = Firework(
[ft,gmx_mdrun_task],
name = 'FORHLR2 GMX mdrun nvt',
spec = {
'_category': 'forhlr2_queue',
'_queueadapter': {
'cpus_per_task': 1,
'ntasks_per_node': 20,
'ntasks': 20,
'queue': 'develop',
'walltime': '00:03:00'
},
'_files_out': {
'log': '*.log',
'trr': '*.trr',
'edr': '*.edr',
'gro': '*.gro' },
'_trackers' : [ gmx_log_tracker ]
} )
pprint(gmx_mdrun_fw.as_dict())
fw_ids = lpad.add_wf(gmx_mdrun_fw)
fw_ids = lpad.get_fw_ids(query={'name':'FORHLR2 GMX mdrun nvt','spec._queueadapter.queue':'develop'})
for fw_id in fw_ids:
try:
print("Deleting {}...".format(fw_id))
lpad.delete_wf(fw_id)
except:
print("Failed deleting {}...".format(fw_id))
ft = FileTransferTask(
mode = 'copy',
files = [ {'src':remote_file, 'dest':os.path.curdir} ] )
gmx_mdrun_task = CmdTask(
cmd = 'gmx',
opt = ['mdrun','-v','-deffnm','nvt'],
stderr_file = 'std.err',
stdout_file = 'std.out',
use_shell = True)
gmx_log_tracker = Tracker('nvt.log')
gmx_mdrun_fw = Firework(
[ft,gmx_mdrun_task],
name = 'FORHLR2 GMX mdrun nvt',
spec = {
'_category': 'forhlr2_queue',
'_queueadapter': {
'cpus_per_task': 1,
'ntasks_per_node': 20,
'ntasks': 40,
'queue': 'normal',
'walltime': '24:00'
},
'_files_out': {
'log': '*.log',
'trr': '*.trr',
'edr': '*.edr',
'gro': '*.gro' },
'_trackers' : [ gmx_log_tracker ]
} )
pprint(gmx_mdrun_fw.as_dict())
fw_ids = lpad.add_wf(gmx_mdrun_fw)
fw_id = list(fw_ids.values())[0]
instead of relying on the returned fw_id, we can also query the Firework added latest
fw_ids = lpad.get_fw_ids(sort=[('created_on',pymongo.DESCENDING)],limit=1)
fw_id = fw_ids[0]
fw_id
We query the remote directory our FireWork ran in
launch_dir = lpad.get_launchdir(fw_id)
launch_dir
c = fabric.Connection('forhlr2') # host defined in ssh config
res = c.run('ls -lht {}'.format(launch_dir)) # look at remote directory contents
glob_pattern = os.path.join(launch_dir,'nvt.*')
res = c.run('ls {}'.format(glob_pattern))
for f in res.stdout.splitlines():
c.get(f)
edr_file = 'nvt.edr'
edr_df = panedr.edr_to_df(edr_file)
edr_df.columns
fig, ax = plt.subplots(3,2,figsize=(10,12))
edr_df.plot('Time','Temperature',ax=ax[0,0])
edr_df.plot('Time','Pressure',ax=ax[0,1])
edr_df.plot('Time','Potential',ax=ax[1,0])
edr_df.plot('Time','Bond',ax=ax[1,1])
#edr_df.plot('Time','Position Rest.',ax=ax[1,1])
#edr_df.plot('Time','COM Pull En.',ax=ax[1,1])
edr_df.plot('Time','Coulomb (SR)',ax=ax[2,0])
edr_df.plot('Time','Coul. recip.',ax=ax[2,1])
mda_trr = mda.Universe('solvated.gro','em_solvated.trr')
# check unique resiude names in system
resnames = np.unique([ r.resname for r in mda_trr.residues ])
resnames
mda_view = nglview.show_mdanalysis(mda_trr)
mda_view.clear_representations()
mda_view.background = 'white'
mda_view.add_representation(repr_type='ball+stick',selection='SDS')
mda_view.add_representation(repr_type='ball+stick',selection='NA')
mda_view.add_representation(repr_type='spacefill',selection='AUM',color='yellow')
mda_view.center()
mda_view
try:
del mda_trr
except:
pass
try:
del mda_view
except:
pass
top = 'sys.top'
gro = 'nvt.gro'
ndx = 'nvt.ndx'
lpad = LaunchPad.auto_load()
gmx_grompp = gromacs.grompp.Popen(
f='npt.mdp',n=ndx,c=gro,r=gro,o='npt.tpr',p=top,
stdout=subprocess.PIPE,stderr=subprocess.PIPE)
out = gmx_grompp.stdout.read()
err = gmx_grompp.stderr.read()
print(err)
print(out)
Utilize fabric to transfer files files to remote resource conveniently:
c = fabric.Connection('forhlr2') # host defined in ssh config
res = c.run('ws_find fw') # get remote directory of Fireworks workspace
res.command
now = datetime.now().isoformat()
remote_path = os.path.sep.join((res.stdout.strip(),'file_transfer',now))
remote_path
res = c.run(' '.join(['mkdir','-p',remote_path]))
file_name = 'npt.tpr'
local_file = os.path.sep.join((prefix,file_name))
remote_file = os.path.sep.join((remote_path,file_name))
res = c.put(local_file,remote_file)
ft = FileTransferTask(
mode = 'copy',
files = [ {'src':remote_file, 'dest':os.path.curdir} ] )
gmx_mdrun_task = CmdTask(
cmd = 'gmx',
opt = ['mdrun','-v','-deffnm','npt'],
stderr_file = 'std.err',
stdout_file = 'std.out',
use_shell = True)
gmx_log_tracker = Tracker('npt.log')
gmx_mdrun_fw = Firework(
[ft,gmx_mdrun_task],
name = 'FORHLR2 GMX mdrun npt',
spec = {
'_category': 'forhlr2_queue',
'_queueadapter': {
'cpus_per_task': 1,
'ntasks_per_node': 20,
'ntasks': 20,
'queue': 'normal',
'walltime': '24:00:00'
},
'_files_out': {
'log': '*.log',
'trr': '*.trr',
'edr': '*.edr',
'gro': '*.gro' },
'_trackers' : [ gmx_log_tracker ]
} )
fw_ids = lpad.add_wf(gmx_mdrun_fw)
fw_id = list(fw_ids.values())[0]
fw_id
instead of relying on the returned fw_id, we can also query the Firework added latest
fw_ids = lpad.get_fw_ids(sort=[('created_on',pymongo.DESCENDING)],limit=1)
fw_id = fw_ids[0]
fw_id
We query the remote directory our FireWork ran in
launch_dir = lpad.get_launchdir(fw_id)
launch_dir
c = fabric.Connection('forhlr2') # host defined in ssh config
res = c.run('ls -lht {}'.format(launch_dir)) # look at remote directory contents
glob_pattern = os.path.join(launch_dir,'npt.*')
res = c.run('ls {}'.format(glob_pattern))
for f in res.stdout.splitlines():
c.get(f)
edr_file = 'npt.edr'
edr_df = panedr.edr_to_df(edr_file)
edr_df.columns
fig, ax = plt.subplots(3,2,figsize=(10,12))
edr_df.plot('Time','Temperature',ax=ax[0,0])
edr_df.plot('Time','Pressure',ax=ax[0,1])
edr_df.plot('Time','Potential',ax=ax[1,0])
edr_df.plot('Time','Bond',ax=ax[1,1])
#edr_df.plot('Time','Position Rest.',ax=ax[1,1])
#edr_df.plot('Time','COM Pull En.',ax=ax[1,1])
edr_df.plot('Time','Coulomb (SR)',ax=ax[2,0])
edr_df.plot('Time','Coul. recip.',ax=ax[2,1])
mda_trr = mda.Universe('nvt.gro','npt.trr')
# check unique resiude names in system
resnames = np.unique([ r.resname for r in mda_trr.residues ])
resnames
mda_view = nglview.show_mdanalysis(mda_trr)
mda_view.clear_representations()
mda_view.background = 'white'
mda_view.add_representation(repr_type='ball+stick',selection='SDS')
mda_view.add_representation(repr_type='ball+stick',selection='NA')
mda_view.add_representation(repr_type='spacefill',selection='AUM',color='yellow')
mda_view.center()
mda_view
substrate = mda_trr.select_atoms('resname AUM')
substrate.masses= ase.data.atomic_masses[ase.data.atomic_numbers['Au']]
substrtate_com_traj = np.array([substrate.center_of_mass() for ts in mda_trr.trajectory ])
substrtate_rgyr_traj = np.array([substrate.radius_of_gyration() for ts in mda_trr.trajectory ])
fig = plt.figure(figsize=(12,10))
ax = fig.add_subplot(111, projection='3d',azim=-30)
ax.plot(*substrtate_com_traj.T)
ax.scatter(*substrtate_com_traj[0,:],color='green')
ax.scatter(*substrtate_com_traj[-1,:],color='red')
plt.plot(substrtate_rgyr_traj)
try:
del mda_trr
except:
pass
try:
del mda_view
except:
pass